DATA624 HW4

(a) Using visualizations, explore the predictor variables to understand their distributions as well as the relationships between predictors.

kable(summary(Glass[,1:5]))%>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

kable(summary(Glass[,6:10]))%>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))

There 9 numerical predictors (RI, Na, Mg, Al, Si, K, Ca, Ba, and Fe) and one categorical target with 6 levels (Type). The data set is highly skewed toward categories 1 and 2 with 146 out of 214 observations falling into those two categories alone. Ba, and Fe have their first quartile value at 0 indicating that they have a lot of zeros in their distributions.

Glass[,1:9] %>%
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram(fill = '#4575b4') +
  theme(panel.background = element_blank(), legend.position="top")

Histograms of the predictors show that all of them except maybe one, Si are moderately to highly skewed. The histogram below has been colored by the target categories and shows that the distributions for different targets values are slightly different (especially for Al and Na) indicating that those variables may be good predictors.

Glass %>%
  gather(-Type, key = "var", value = "val") %>%
  ggplot(aes(x = val, fill=Type)) +
  geom_histogram(bins=10, alpha=1) +
  facet_wrap(~ var, scales = "free") +
  scale_fill_manual("target",
                    values = c('#d73027','#fc8d59','#fee090','#e0f3f8','#91bfdb','#4575b4')) +
  xlab("") +
  ylab("") +
  theme(panel.background = element_blank(), legend.position="top")

Those differences in the distributions of the predictors when separated by the target variable are even more evident in the bar plots below. All of the predictors show a lot of variation in their distributions for each target value. We can also see a lot of outliers which are indicated by the red dots to the top and bottom of the plot ‘whiskers’.

Glass %>%
  gather(-Type,key = "var", value = "val") %>%
  ggplot(aes(x=factor(Type), y=val)) +
  geom_boxplot(width=.5, fill="#58BFFF", outlier.colour="red", outlier.size = 1) +
  stat_summary(aes(colour="mean"), fun.y=mean, geom="point",
               size=2, show.legend=TRUE) +
  stat_summary(aes(colour="median"), fun.y=median, geom="point",
               size=2, show.legend=TRUE) +
  facet_wrap(~ var, scales = "free", ncol=3) +
  labs(colour="Statistics", x="", y="") +
  scale_colour_manual(values=c("#9900FF", "#3300FF")) +
  theme(panel.background=element_blank())

Only RI and Ca seem to be very highly correlated with each other with a correlation of about 0.81.

corrs <- round(cor(Glass[,1:9]),3)
ggcorrplot::ggcorrplot(corrs,
       type = 'lower', lab=T, lab_size=2)

Glass %>% 
  ggscatmat(color="Type", alpha=0.6) +
  scale_color_manual(values=c('#d73027','#fc8d59','#fee090','#e0f3f8','#91bfdb','#4575b4')) +
  theme(panel.background=element_blank(), legend.position="top",
        axis.text.x = element_text(angle=-40, vjust=1, hjust=0))

highCorr <- findCorrelation(corrs, .75)
highCorr

## [1] 7

The findCorrelation function suggests removing column 7 which corresponds to Ca.

There are no missing values in the dataset.

missing <- data.frame(t(apply(is.na(Glass), 2, sum)))
# kable(missing, col.names = "NA's")
row.names(missing) <- c("NA's")
kable(missing) %>%
  kable_styling(bootstrap_options = c("condensed"))

(b) Do there appear to be any outliers in the data? Are any predictors skewed?

As noted above most if not all of the predictors have both outliers and skewed distributions as well as a large percentage of zeroes in a few predictors.

(c) Are there any relevant transformations of one or more predictors that might improve the classification model?

skewness_statistic <- apply(Glass[,1:9], 2, skewness)
ratio_max_to_min <- apply(Glass[,1:9], 2, max)/apply(Glass[,1:9], 2, min)
kable(data.frame(skewness_statistic, ratio_max_to_min)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Because there are so many zero values in the data it’s hard to estimate the severity of the skew using the max to min ratio since division by zero is undefined. But we clearly have some large skewness values and can clearly see the skew in the histograms.

The skewed distributions might be improved by log, square root, or inverse transformations or a Box-Cox transformation.

Unfortunately because of the zero values Mg, K, Ba, and Fe cannot be transformed using Box-Cox. I’m still not sure how to handle that although I found this post by Rob Hyndman on the subject… Transforming data with zeros.

# https://stackoverflow.com/questions/46485024/how-to-use-boxcoxtrans-function-in-r
GlassBC <- apply(Glass[,1:9], 2, BoxCoxTrans)
trans_Glass <- purrr::map2(GlassBC, Glass[,1:9], function(x, y) predict(x, y)) 
trans_Glass_df = as.data.frame(do.call(cbind, trans_Glass))
kable(head(trans_Glass_df)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

trans_Glass_df %>%
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram(fill = '#4575b4') +
  theme(panel.background = element_blank(), legend.position="top")

I’m not seeing a lot of improvement in the histograms after transformation with Box-Cox so let’s try something else… Another post suggests three ways of handling the zeros:

1. Add a constant value © to each value of variable then take a log transformation
   
2. Impute zero value with mean
   
3. Take square root instead of log for transformation²

Let’s try the square root transformation first.

# Square Root transformation
Glass_SqRt <- data.frame(apply(Glass[,1:9], 2, sqrt))

Glass_SqRt %>%
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram(fill = '#4575b4') +
  theme(panel.background = element_blank(), legend.position="top")

The square root transformation seems to have done a little bit better at taming our distributions, however they are still highly skewed in some cases.

Let’s try one more of the suggestions by using \(log(x+1)\) to transform our data.

# log(x+1) transformation
Glass_log <- data.frame(apply(Glass[,1:9]+1, 2, log))

Glass_log %>%
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram(fill = '#4575b4') +
  theme(panel.background = element_blank(), legend.position="top")

Maybe a little better, but still not great!

kable(nearZeroVar(Glass, names = TRUE, saveMetrics=TRUE)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

The caret function nearZeroVar does not indicate that the variabes with high frequenciees of zeroes should be removed so another solution might be to try the alternative Box-Cox transformation that Hyndman suggested in his blog post.³

Alt_BoxCox_Glass <- log1p(Glass[,1:9])
kable(head(Alt_BoxCox_Glass))  %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

Alt_BoxCox_Glass %>%
  gather() %>% 
  ggplot(aes(value)) +
    facet_wrap(~ key, scales = "free") +
    geom_histogram(fill = '#4575b4') +
  theme(panel.background = element_blank(), legend.position="top")

I think that in this case it mght be necessary to know if these are real zeroes, or ‘censored’ data where the zeroes may actually be trace amounts that are below the detection limit. In this case Kuhn and Johnson state:

when a sample has a value below the limit of detection, the actual limit can be used in place of the real value. For this situation, it is also common to use a random number between zero and the limit of detection.⁴

(a) Investigate the frequency distributions for the categorical predictors. Are any of the distributions degenerate in the ways discussed earlier in this chapter?

# reorder columns by number of levels in factors
Soybeans <- Soybean[, c(1,2,7,8,22,23,29,30,4,5,9:11,14:16,19,25,27,36,
                       3,6,12,13,17,18,20,21,24,26,28,31:35)]
kable(summary(data.frame(Soybeans[,1:2]))) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

kable(summary(Soybeans[,3:8])) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

kable(summary(Soybeans[,9:14])) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

kable(summary(Soybeans[,15:20])) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

kable(summary(Soybeans[,21:28])) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

kable(summary(Soybeans[,29:36])) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

par(mar=c(1.2,1.2,1.2,1.2))
par(mfrow=c(9,4))
for (i in 1:ncol(Soybeans)) {
  barplot(summary(Soybeans[,i]), main=names(Soybeans[i]))}

(b) Roughly 18 % of the data are missing. Are there particular predictors that are more likely to be missing? Is the pattern of missing data related to the classes?

plot_missing(Soybeans) 

gg_miss_var(x = Soybeans, facet = Class)

missingness <- Soybeans %>%
  group_by(Class) %>%
  miss_var_summary()
missingness

nearZeroVar(Soybeans)

## [1] 17 30 31

kable(nearZeroVar(Soybeans, names = TRUE, saveMetrics=TRUE)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

(c) Develop a strategy for handling missing data, either by eliminating predictors or imputation.